Compute first guess

ABSTRACT

THE PRESENT INVENTION GENERALLY RELATES TO A NOVEL METHOD FOR OPERATING AN ELECTRONIC COMPUTER TO PERFORM INTEGRATION OPERATIONS UPON A SET OF ORDINARY DIFFERENTIAL EQUATIONS. THE OVERAL METHOD IS ILLUSTRATED IN THE FLOW CHART OF FIG. 1. ACCORDING TO A FIRST ASPECT OF THE INVENTION OCCURRING IN BLOCK 7 OF FIG. 1, THERE IS DISCLOSED AN AVERAGING PROCESS COMPRISING THE FOLLOWING STEPS. SELECT AN ORDER OF ACCURACY P FOR THE UNDERLYING LINEAR MULTISTEP FORMULAE SOLUTIONS WHICH ARE TO BE AVERAGED. SELECT A NUMBER D BY WHICH THE AVERAGING PROCESS IS TO RAISE THE ORDER OF ACCURACY P OF THE INDIVIDUAL UNDERLYING SOLUTIONS. SELECT A LINEAR MULTISTEP FORMULA OF ORDER P AND A SUFFICIENTLY LARGE STEP NUMBER K TO PRODUCE UNDERLYING SOLUTIONS, WHEREIN THE FORMULA MUST CONTAIN   AT LEAST M$D FREE PRIMARY PARAMETERS AND ONE SECONDARY PARAMETER C, THE COEFFICIENT OF THE LEADING DERIVATIVE FOUND IN THE FORMULA. DETERMINE THE NUMBER N AND FORM OF THE CONSTRAINTS WHICH THE WEIGHT MUST SATISFY FOR THE AVERAGING PROCESS TO RAISE THE ORDER FROM P TO P+D. SELECT N DIFFERENT SETS OF VALUES OF THE FREE PARAMETERS. COMPUTER THE N WEIGHTS SATISFYING THE N CONSTRAINTS. COMPUTE N SOLUTION OF THE DIFFERENTIAL EQUATIONS, EACH OF ORDER P, EACH ASSOCIATED WITH ONE OF THE PARAMETER SETS SELECTED, AND AVERAGE THESE N SOLUTIONS TO OBTAIN A SOLUTION HAVING AN ORDER OF ACCURACY OF P+D. ACCORDING TO A SECOND ASPECT OF THE INVENTION WHICH OCCURS IN BLOCK 3 OF FIG. 1, A-STABILITY IS MAINTAINED BY CARRYING OUT THE FOLLOWING STEPS. FIRST, FIND THE FORM OF THE A-STABILITY CONSTRAINTS TO BE IMPOSED ON THE PARAMETERS OF THE FORMULA SELECTED IN THE AVERAGING PROCESS. THESE CONTRAINTS CONSIST OF INFINITELY MANY LINEAR, AND FINITELY MANY NONLINEAR INEQUALITIES. NEXT, FIND ALL SETS OF PARAMETER VALUES SATISFYING ALL LINEAR A-STABILITY CONSTRAINTS. FROM THESE SETS OF PARAMETER VALUES SELECT, BY METHODS OF NONLINEAR PROGRAMMING, SET SATISFYING THE NONLINEAR CONSTRAINTS AS WELL. IT SHOULD BE NOTED THAT THE ABOVE SELECTION OF PARAMETERS IS CONTINGENT UPON SELECTING AN APPROPRIATE LINER MULTISTEP FORMULA AS SET FORTH IN THE FIRST ASPECT OF THE INVENTION. IN OTHER WORDS, IF A PROPER LINER MULTISTEP FORMULA MEETING THE SPECIFIED CONSTRAINTS IS CHOSEN, IT IS POSSIBLE TO ALSO SELECT PROPER PARAMETER VALUES TO OBTAIN THE DESIRED STABILITY RESULTS THE THIRD AND MOST SIGNIFICENT ASPECT OF THE PRESENT INVENTION COMPRISES THE ACTUAL INTEGRATION PROCEDURE. THIS OPERATION OCCURS PRIMARILY IN BLOCK 7 OF FIG. 1. AN INITILIZATION STEP INCLUDES LOADING THE SYSTEM WITH INFORMATION PERTAINING TO THE PREVIOUSLY SELECTED LINEAR MULTISTEP FORMULA AS WELL AS WITH THE SOLUTION VALUES NECESSARY FOR STARTING THE INTEGRATION. ALSO, THE INTEGRATION STEP H, STEP NUMBER K AND THE PREVIOUSLY ENUMERATED PARAMETER VALUE SELECTIONS ARE LOADED DURING INITILIZATION. THE FIRST STEP OF THE ACTUAL INTEGRATION PROCESS IS TO UPDATE THE COUNTER N OF INTEGRATION STEPS, A TEST IS THEN CARRIED OUT TO SEE IF ALL REQUIRED STEPS HAVE BEEN COMPLETED. IF YES, THE AVERAGED SOLUTION Z IS PRINTED OUT AND THE PROGRAM ENDS. IF NOT, THEN AN INTEGRATION STEP IS PERFORMED, COMPRISING COMPUTING T HE FIRST UNDERLYING SOLUTION X(1) OF THE SYSTEM OF ORDINARY DIFFERENT EQUATIONS BY THE NEWTON-RAPHSON METHOD, AS WELL AS AUXILIARY QUA

v UNITED STATES PATENT OFFICE Published at the request of the applicant or owner in accordance with the Notice of Dec. 16, 1969, 869 0.G. 687. The abstracts of Defensive Publication applications are identified by distinctly numbered series and are arranged chronologically.

The heading of each abstract indicates the number of pages of specification, including claims and sheets of drawings contained in the application as originally flied. The files of these applications are available to the public for inspection and reproduction may be purchased for 30 cents a sheet.

Defensive Publication applications have not been examined as to the merits of alleged invention. The Patent Otiice makes no assertion as to the novelty of the disclosed subject matter.

PUBLISHED JULY 24, 1973 T912 013 HIGHLY ACCURATE UNCONDITIONALLY E G STABLE METHOD FOR INTEGRATING 1 ORDINARY DIFFERENTIAL EQUATIONS Werner Liniger and Farouk M. Odeh, Yorktown Heights, p 0F N N .Y., assigiiors :3 Ir11enational Business Machines Cor- M U Ln ST EP so LUT l 0 NS WH CH oration, rmo p Filed June 30, 1972, Ser. No. 267,806 ARE TO BE AVERAGE D Int. Cl. G06f 15/32 2 US. Cl. 4441 8 Sheets Drawing. 52 Pages Specification SELECT N B ER d Y w H I c H The present invention generally relates to a novel method for operating an electronic computer to perform IS TO SE THE integration operations upon a set of ordinary differenp NG so T tial equations. The overall method is illustrated in the flow chart of FIG. 1. According to a first aspect of the 3 t invention occurring in block 7 of FIG. 1, there is disclosed an averaging process comprising the following FORMULA 0F p: steps. Select an order of accuracy p for the underlying k [0 linear multistep formulae solutions which are to be averaged. Select a number d by which the averaging process is to raise the order of accuracy p of the individual underlying solutions. Select a linear multistep formula of order p and a sufiiciently large step number k to produce 6, underlying solutions, wherein the formula must contain at least mzd free primary parameters and one secondary parameter c, the coeflicient of the leading derivative TERM IN TH E FORMULA found in the formula. Determine the number N and form 4 of the constraints which the weights must satisfy for the averaging process to raise the order from p to p+d. Se-

lect N difierent sets of values of the free parameters. WH|QN WEIGHTS Compute the N weights satisfying the N constraints. Compute N solutions of the differential equations, each of or- F l der p, each associated with one of the parameter sets A selected, and average these N solutions to obtain a so- RMSE FROM p To p+d lution having an order of accuracy of p+d. 5 4' According to a second aspect of the invention which occurs in block 3 of FIG. 1, A-stability is maintained by N DI SETS 0F carrying out the following steps. First, find the form of the A-stability constraints to be imposed on the parameters of the formula selected in the averaging process.

These constraints consist of infinitely many linear, and 6 finitely many nonlinear inequalities. Next, find all sets of COM PUTE N WEIGHTS BY SATISFY'NG parameter values satisf in all linear A-stability conin s, Y g N CONSTRAINTS BLOCK 4 (THESE From these sets of parameter values select, by meth- REPRESENT N L|NEAR ALGEBRMC Ode of nonlinear programming, sets satisfying the nonlinear constraints as well. It should be noted that the above selection of parameters is contingent upon selecting an appropriate linear multistep formula as set forth in the first aspect of the invention. In other words, if a N proper linear multistep formula meeting the specified constraints is chosen, it is possible to also select proper EQUATIONS! EACH ORDER- EACH parameter values to obtain the desired stability results. AT'ED WITH ONE PARAMETER The third and most significant aspect of the present invention comprises the actual integration procedure. This N L U 0N3 operation occurs primarily in block 7 of FIG. 1. An initialization step includes loading the system with infor- To OBTAIN SOLUTION 0F ORDER p PUBLISHED J ULY Ql, 1973 mation pertaining to the previously selected linear multistep formula as Well as with the solution values necessary for starting the integration. Also, the integration step h, step number k and the previously enumerated parameter value selections are loaded during initialization. The first step of the actual integration process is to update the counter 21 of integration steps, a test is then carried out to see if all required steps have been completed. If yes, the averaged solution z is printed out and the program ends. If not, then an integration step is performed, comprising computing the first underlying solution x of the system of ordinary differential equations by the NeWton-Raphson method, as well as auxiliary quantities associated with said first underlying solution.

W. LINIGER ET July 24, 1973 HIGHLY ACCURATE UNCONDITIONALLJY STABLE METHOD FOR TNTEGRATING ORDINARY DIFFERENTIAL EQUATIONS Filed June 30, 1972 FIG. I

SELECT NUMBER d BY VIHICH AVERAGING IS TO RAISE THE ORDER p OF THE UNDERLYING SOLUTIONS 8 Sheets-Sheet 1 FIG. 2

A-STABILITY SELECT FORMULA OF ORDER p, SUFFICIENTLY LARGE STEP NUMBER k,TO PRODUCE UNDER- LYING SOLUTIONS FORMULA MUST CONTAIN AT LEAST d FREE PRIMARY PARAMETERS AND ONE SECONDARY PARAMETER 0, THE

COEFFICIENT OF THE LEADING DERIVATIVE TERM IN THE FORMULA DETERMINE NUMBER N AND FORM OF CONSTRAINTS YIHICH VIEIGHTS MUST SATISFY FOR AVERAGING TO RAISE ORDER FROM p TO p+d SELECT u DIFFERENT SETS oF VALUES 0F THEm FREE PARAMETERS COMPUTE N VIEISHTS BY SATISFYING N CONSTRAINTS BLOCIII (THESE REPRESENT N LINEAR ALGEBRAIC EQUATIONS FOR VIEIGHTS COMPUTE N SOLUTIONS OF DIFFERENTIAL EQUATIONS, EACH OF ORDER p,EACH

ASSOCIATED YIITH ONE PARAMETER SET. AVERAGE THESE N SOLUTIONS TO OBTAIN SOLUTION OF ORDER p+d.

FIND ALL SETS OF PARAMETER VALUES SATISFYING ALL LINEAR ASTABILITY CONSTRAINTS. (SELECTION OF FORMULA IN FIG.I BLOCK 5 AND OF VALUE OF SECONDARY PARAMETERS MUST GUARANTEE EXISTENCE OF SUCH SETS) \FROM THESE SETS OF PARAMETER VALUES,

EXISTENCE OF SUCH SETS.)

THE N SETS OF PARAMETER VALUES NEEDED IN AVERAGING PROCESS ARE ALL CHOSEN FROM SETS DEFINED IN BLOCKS THEN ALL UNDERLYING SOLUTIONS WILL BE A-STABLE.

AVERAGE A-STABLE UNDERLYING SOLUTIONS TO PRODUCE A-STABLE AVERAGED SOLUTION OF HIGHER ORDER July 24, 1973 w. LINIGER ET AL T912013 HIGHLY ACCUIIA'TTJ' UNCONDITIUNALLY STATSTITLI MIEI'I'ITOI) T'OII IN'IFIITIIA'T'TNI'I ORDINARY DIFFERENTIAL EQUATIONS Filed June 30, 1972 8 Sheets-Sheet 2 (FIGI BLOCK T) INITIALIZATION UPDATE COUNTER 2 OF INTEGRATION STEPS n YES I PRIIITIIUT 5 m AVERAGED N UPDATE UNDERLYING SOLUTION x ,USING SOLUTION 1 IIEIIIToII-RAPIIsoII IIETIIoII. UPIIATE AUXILIARY oUAIITITIEs ASSOCIATED IIIITII x END UPDATE AVERAGED SOLUTION z UPDATE PERTURBATION BEWEEN SOLUTION I I AND XU, BY LINEARIZATION AROUND I UPDATE AUXILIARY oUAIITITIEs ASSOCIATED IIIITIIE y 1973 w. LINIGER ET T911013 HIGHLY ACCURATE UNCONDITIONALLY STABLE METHOD FOR INTEGRATING ORDINARY DIFFERENTIAL EQUATIONS Filed June 30, 1972 8 Sheets-Sheet 5 FIG. 4 NONLINEAR PROCEDURE 5 FOR UPDATING FIRST UNDERLYING SOLUTION LINEAR PROCEDURE FOR UPDATING UNDERLYING I FIG BLOCK 5) SOLUTIONS x (P) ,P= 2---,N.

COIIPUTE FIRST GUESS 'ITHHBY POLYNOMIAL EXTRAPOLATION 0F SUFFICIENTLY HIGH DEGREE COIIPUTE F J (FIG?) BLOCK 8 I CONPUTE Fl RST GUESS OF (p) BY 2 PERTURBATION POLYNOM I AL EXTRAPOLATION SOLVE FOR I FIRST AND ONLY NEWTON -RAPHSON COR R ECT I ON (I) 9 (I I I IREQU'RES ONE CONPUTE NEWTON-RAPHSON L/U FACTORIZAT'OIII 2 CORRECTION OF PERTURBATION BACK SUBSTIIUTIONS A ASSOCIATED WITH X) a f COIIPUIE FINAL UPDATED VALUE v OF FIRST UNDERLYING SOLUTION UPDATE PERTURBATIQN (I)+=,(I)+ (1)+ gIPH EIPH gIPI+ RECOIIPUTE f AUXILIARY QUANTITIES July 24, 1973 HIGHLY ACCURATE UNCONDITIONALLY STABLE METHOD FOR INTEGRATING Filed Junb 50, 1972 w. LINIGER ET AL ORDINARY DIFFERENTIAL EQUATIONS FIG. 6

AVERAGING (FIGS BLOCK IO) 8 Sheets-Sheet 4 UPDATE AVERAGED SOLUTION USING UPDATED SOLUTION X(1)+AND UPDATED FIG. 7

P=2 AVERAGING PROCESS 3 d=2 I 82 21, T5 T, $5 T GOMPUTE 5 SOLUTIONS OF ORDER 2,ASSOG|ATEO WITH (r,s)=(T,2),(5,2),AND(7,1) RESPEGTIVELY. AVERAGE 3 SOLUTIONS TO OBTAIN SOLUTION OF ORDER 4 July 24, 1973 w. LINIGER ET AL T912,0l3

HIGHLY ACCURATE UNCONDITIONALLY STABLE METHOD FOR INTEGRATING ORDINARY DIFFERENTIAL EQUATIONS Filed June 30, 1972 8 Sheets-Sheet 5 FIG. 8

A-STABI IT Y (SAME AS BLOCK I FIC 2 I FOR THE VALUE OF SECONDARY PARAMETER CHOSEN IN FIC.T

2 BLOCK3 THESE SETS FORM. THE ELLIPTICAL REGION OF FICIS PLUS REGION INSIDE TANCENTS I AND I THE NONLINEAR CONSTRAINT DEFINES REGION ABOVE HYPERBOLA PASSING THROUGH P2,P4,P6,P7 AND P1.ALL SETS 0F PARAMETER VALUES SATISFYING CONSTRAINTS 0F BLOCK 2 ABOVE ALSO STATISFY THIS CONSTRAINT SETS DEFINED IN BLOCKS 0F FIC.T 4V BELONC TO RECIONS DEFINED IN BLOCKS 2AND 3 ABOVE 5 (SAME AS BLOCK 5 FIG. 7)

July 24, 1973 w. LINIGER E T912913 HIGHLY ACCURATE UNCONDITIONALLY STABLE METHOD FOR INTEGRATING ORDINARY DIFFERENTIAL EQUATIONS Filed June 30, 1972 8 Sheets-Sheet 6 F IG. 9

INTEGRATION PROCEDURE INITIALIZATION [11-[52] NT UPDATE STEP COUNTER 2 ALL STEPS CARRIED OUT? "PRINT OUT" 5/ UPDATE SOLUTION #1 AVERAGED SOLUTIONI [55-sT,Us-aT] 6\ lNITlALlZE SOLUTION COUNTER n 1 [To] END 10 UPDATE ALL (5) AVERAGED SOLUTIONS SOLUTION 1-. UPDATED P 5? UPDATE SOLUTION /9 COUNTER p: [D1] Filed June 30, 1972 July 24, 1973 w. LINIGER ETAL T912,013

HIGHLY ACCURATE UNCONDITIONALLIY STABLE METHOD FOR INTEGRATING ORDINARY DIFFERENTIAL EQUATIONS 8 Sheets-Sheet '7 FIG. I0 FIG. 11

NONLINEAR PROCEDURE LINEAR PROCEDURE FOR UPDATING FIRST FOR UPDATING UNDERLYING SOLUTION UNDERLYING SOLUTIONS X 2,3

COMPUTE FIRST GUESS T [56,51] 3"[58] if? [59] COMPUIE FIRST GUESS FOR PERTURBATION g SOLVE FOR NEIIITON- RAPHSON conmnonefl; [so-e5] (L\U FACIORIZATIOII 1 [eI],IIII0 IIAcII SU-BSTITUIIONS 1 [e2 ,63]

COMPUTE NEWTON-RAPHSON CORRECTIONS T0 PERTURBATION 6(P)+ [Is 14] BY I BACK SUBSIITUTIONSI [mu] COMPUTE FINAL UPDATED VALUE OF FIRST soIuII0II:I' 4] UPDATE PERTURBATION AUXILIARY QUANTITIES ITS-8Q] RECOMPUIE I AUXILIARY ouAIIIIIIEs: [TB-87] FIG. I2

AVERAGING UPDATE AVERAGED SOLUTION! E38] (SAME As FIG.6 I

July 24, 1973 w. LINIGER ET AL T912013 HIGHLY ACCURATE UNCONDITIONALLY STABLE METHOD FOR INTEGRATING ORDINARY DIFFERENTIAL EQUATIONS Filed June 30, 1972 8 Sheets-Sheet 8 Fl G. E

ASYMPTOTE 1 

